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THE  VISCOELASTIC  DEFLECTION 
OF  AN  INFINITE  FLOATING  ICE  PLATE 
SUBJECTED  TO  A CIRCULAR  LOAD 

Shunsuke  Takagi 


INTRODUCTION 

Since  ancient  times  floating  ice  plates  have  been  used  to  cross  rivers  and  lakes.  During  recent 
years  traffic  load  on  frozen  rivers  and  lakes  has  greatly  increased,  and  at  the  same  time  vehicles  have 
become  heavier.  Aircraft  landing  and  parking  facilities  also  have  added  loads  on  these  bodies  of  water. 
In  addition,  during  the  past  several  years,  oil  companies  have  started  to  use  ice  plates  as  drilling  plat- 
forms. Thus,  we  now  need  to  acquire  a more  detailed  understanding  of  the  creep  of  ice  plates. 

Formulation  of  the  creep  of  a floating  ice  plate  began 
after  World  War  II  with  the  intense  development  of  the  linear 
viscoelasticity  theory.  In  1947  Golushkevich  (referred  to  by 
Kheysin'®)  presented  an  analysis  assuming  that  ice  behaves 
elastically  for  volumetric  deformations  and  viscoelastically 
for  deviatoric  deformations.  Kheysin*®  used  a general  visco- 
elastic thin-plate  theory  to  analyze  the  infinite  floating  ice 
plate.  He  used  the  Maxwell  unit  (Fig.  1)  only,  and  considered 
only  a concentrated  load.  Nevel"  also  used  the  Maxwell  unit 
only,  but  considered  a distributed  load.  He  limited  his 
numerical  computation  only  to  the  center  of  the  load. 

William  L.  Ko,  as  reported  by  Garbaccio,'*  * used  the 
Maxwell-Voigt  type  four-element  model  (Fig.  1),  which  is 
known  to  represent  the  creep  of  ice  satisfactorily  (Jellinek 
and  Brill*).  In  addition  to  thin-piate  theory,  Ko  used 
Reissner’s  plate  theory,  which  includes  the  deflection  due 
to  vertical  shear  forces.  Garbaccio®  numerically  evaluated 
Ko’s  solution  for  specific  values  of  material  constants  rather 
than  for  nondimensional  parameters.  Garbaccio’s  numerical  answers  show  that  the  discontinuity  of 
the  load  distribution  yields  a strong  influence  on  the  values  of  deflection.  It  is  reasonable  to  suspect 
that  his  numerical  evaluation  may  contain  some  errors. 

lAkunin*  has  solved  the  same  problem  as  Ko,  but  he  used  only  thin-plate  theory.  Unfortunately, 
only  an  abstract  of  lAkunin’s  work  is  available  to  western  researchers. 

Katona’  and  Vaudrey  and  Katona'^  solved  the  same  problem  with  a finite-element  viscoelastic 
computer  program. 

We  solved  this  problem  analytically  by  use  of  thin  plate  theory,  and  also  developed  an  effective 
method  of  numerical  integration  of  the  solution  integrals.  However,  the  theoretical  curves  did  not 
satisfactorily  fit  the  field-test  curves.  It  is  now  evident  that  a large  scale  laboratory  test  eliminating 
the  variation  due  to  natural  conditions  must  be  carried  out  and  the  theoretical  assumptions  must  be 
tested. 


Moxwel  I 
Unit 


E2  Voigt 
Unit 


Figure  1.  Maxwell-Voigt  type 
four  element  model. 


1.  THE  PROBLEM 


We  shall  consider  the  viscoelastic  ice  plate  floating  on  water  extending  horizontally  to  infinity. 

We  shall  use  the  Maxwell-Voigt  type  four-element  model  (Fig.  1 ) to  describe  the  viscoelastic  deforma- 
tion of  ice. 

Using  the  notation  of  Fig.  1,  we  can  show  that  this  model  gives  the  stress-strain  relationship 
which  we  show  in  an  operator  form, 


e = 


(1.1) 


where  t is  time.  To  extend  the  one-dimensional  relationship  (1 .1 ) to  the  three-dimensional  relation- 
ship, we  assume,  as  explained  by  Fliigge,*  that  e and  a are  deviatoric  and  relate  them  by 


a = 2Ge 

where  G is  the  rigidity  modulus  relative  to  the  three-dimensional  deformation.  Using  (1.1),  2G  is 
given  as  an  operator 


J_ 

2G 


_L 


(1.2) 


The  differential  equation  describing  the  deflection  w of  an  elastic  plate  on  water  is 


Dy^w  + pw  = q (1 .3) 

where  is  the  biharmonic  operator 


p the  density  of  water,  q the  load  per  unit  area,  D the  flexural  rigidity  defined  by 

D = 2Gh^|[^2(^ -v)]  (1.5) 

in  which  h is  the  thickness  of  the  ice  plate,  and  v Poisson’s  ratio.  Substituting  2G  from  (1 .2)  into 
(1.5),  and  D thus  found  into  (1.3),  we  find  the  differential  equation  governing  the  viscoelastic  de- 
flection of  a floating  ice  plate.  We  shall  show  this  equation  later  in  the  nondimensional  form. 

We  assume  the  load  q to  be  a step  loading  applied  at  r = 0 and  distributed  uniformly  over  a circle 
of  radius  a with  the  center  at  origin.  Then,  letting  r be  the  radial  distance  from  origin 

q = qo  U(t)  for  0 S r < (7 

(1.6) 

= 0 for  <j  < r 

where  U(t)  is  the  step  function,  and  t the  time.  Our  problem  is  axisymmetric,  and  the  biharmonic 
operator  V*  reduces  to 
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We  shall  nondimensionalize  our  differential  equation.  We  define  the  characteristic  length  8 by 


S*  = Eoh^l[^2fi0  -•')] 


J-  = J-  + J- 


We  have  chosen  Eq,  rather  than  E^  or  £’21  to  define  8,  because  Eq  is  related  to  the  secondary 
creep  (Nevel**),  which  is  the  main  interest  in  our  field  observation. 

Let  D,  be  defined  by 

D,  = (1.S 

Use  of  (1 .4)  and  (1 .7)  changes  (1 .9)  to 

O]  = 2G/Eq.  (1.1 

Substituting  G in  (1.2),  (1.10)  becomes 


o,  = 1/  / + -V  + ^ 


We  choose  nondimensional  time  T 


T = Eq  r/rj, 
and  a parameter  r 

T = t?!  f 2/(^2  ^o)- 
Then  (1.11)  becomes 

' dT  bT 


E = Eq/E,  . 


It  is  noted  that 


0 < £ < 1 . 


Clearing  the  denominator,  (1.14)  becomes 


where  use  is  made  of  the  relation 


ET+T?l/t?2  = T 

which  can  be  proved  by  use  of  (1.13),  (1.15),  and  (1.8). 

We  define  the  nondimensional  length  /?  by 

R = r/e.  (1.18) 

We  replace  D in  (1 .3)  with  D in  (1 .9),  and  (1 .3)  becomes 

D■^\7^w+w  = q/p  (1.19) 

where 

With  £>,  given  by  (1.17),  (1.19)  is  the  differential  equation  to  be  solved. 


2.  THE  SOLUTION 

We  denote  the  Hankel  transform  of /"(/?)  by  f(p) 

oo 

m = f f(R)jomRdR 
0 

and  the  two-sided  Laplace  transform  (Van  der  Pol  and  Bremmer'*)  offf(7)  by  g{s) 

^(s)  = 5 J 9(T)e-STdr. 

Wc  denote  the  inverse  of  (2.2)  by 
g(T)  = L-'  [9(5)1. 

Applying  these  two  transforms,  (1.19)  becomes 
D,  vv  + w = q/p 

where 

D,  = . 

fS2  + (1  + t)S  + T 

Applying  the  two  transforms  to  q defined  by  (1 .6),  we  get 
(1/P)t  = [PI{irApi^)](im 

where 


(2.1) 


(2.2) 


(2.3) 


(2.4) 


(2.5) 
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B 


P = ua^q 


and 


A = a/e. 

Thus  the  transformed  solution  is  given  by 

^ = 4r— r 

nApe^  00 +Djp4) 

Performing  the  Hankel  inverse,  we  find 


^ = “7^  / 7 I 4 
nApe^  ^ 1 + D^0'^ 

Performing  the  Laplace  inverse,  we  find 


To  find /.■'  (1/(1  + /?“’)],  we  compute  the  partial  fraction 

1 1 ^ 1 ^ g"*  (r-a2)  _i ^(t  - «! ) 1 

1 + O,  /}4  5 y/D£SC  S + a2  ~ y/DE^  5 + a, 


where  - a,  and  - 02  are  the  roots  of  the  quadratic  equation 
(£  + 0'*)  52  + (^^4  + 1 + ^)5  + ^ = 0. 

They  are  given  by 

“1  = T/r»  + 1 -f  r T JDESC 

“2  2(^4  + 

where 

DESC  = (T0*  + ^ + t)2  - 47(134  + £) 
which  transforms  to 


(2.6) 


U.7) 


(2.8) 


(2.9) 


(2.10) 


(2.11) 


(2.12) 


= [t(^  + 1)-1]2+4t(1-E),  (2,13) 

From  (2.1 3),  it  is  clear  that 

DESC  > 0.  (2  14) 

The  roots  a,  and  a2  are  therefore  always  real.  Moreover,  inspection  of  (2.1 1)  and  (2.12)  shows 
that  both  a,  and  a2  are  always  positive.  Thus  we  find  that 


5 


r 


L-U  ' ] . 

s/DESC  s/DESC 

Substituting  (2.15)  into  (2.9),  the  solution  for  w is  found: 

M 

/ 


(2.15) 


^^^(7^2)  .-g^T  ,.-a,T 

y/DESC  y/DESC 


J,m  dp.  (2.16) 


itApi^  0 
The  radial  and  hoop  stresses  are  given  by 
_ 6D  ( d^w  ^ V dw\ 

^ ' /»2  \ 3r2  ^ 3/-  y 

„ =_60/l|vv 
® a/-2  / 


respectively.  Changing  D to  D,  by  use  of  (1 .9)  and  r to  nondimensional  R by  use  of  (1 .1 8),  they 
become 


, _ 6o£2  n /b'^w  3w  \ 

a - 6pe2  u (]_  3vv  3ijy\ 


where  £>,  is  the  operator  on  T given  by  (1 .1 7).  The  two-sided  Laplace  transform  yields 


. 6o82  / 02  ^ 


® /»2 


2 \r  dR  sRlJ 


(2.17) 

(2.18) 


where  w is  the  Laplace  transform  of  w. 
Using  (2.8)  one  gets 


D■^w  = 


itApi'^  0 
The  Laplace  inverse  of  this  is 


DfW  = 


7r/lp82  Q 

To  find  L"’  [Di/(1  + |3^)] , we  compute  the  partial  fraction, 


X 


1 ^ ’■-“l  _1  ^•-«2  1 

5 1 + zy,  ^ ~ \/DEsc  5 + tti  ~ vsnc  s + <*2  ■ 
Thus  we  find 

- ’•-“1  ,-^.T  ^-«2 

Vl+D,^/  \/5£SC  y/DESC 


Thus  the  inverse  of  (2.1 7)  is 


a,  = / J^ifiA) 


' irAh^ 


JoW- 


1 -V 

0R 


VS£5C 


The  inverse  of  (2.1 8)  is 


6/> 

■nAh^ 


f J^ifiA) 


(^/?) 


(r-«,)e-“l^-(r-a2)e-°2T 


s/DESC 


dfi.  (2.20) 


Tabulation  of  and  Og  becomes  easier  if  linear  combinations  of  (2.1 9)  and  (2.20)  that  do  not  con- 
tain i>  are  computed. 


3.  METHOD  OF  NUMERICAL  INTEGRATION 

It  is  impossible  to  analytically  integrate  the  solution  integrals  (2.16),  (2.19)  and  (2.20).  (Sec 
App.  I.) 

The  direct  numerical  integration  is  inconvenient  because  of  the  slow  convergence  of  the  Bessel 
functions  for  large  values  of  the  independent  variable  We  shall  choose  finite  ranges  of  integration 
that  give  sufficiently  close  approximations.  The  essence  of  our  method  consists  of  the  following 
integration  procedure: 

Consider  the  integral 

/=  f m hm  hm  o.i) 

0 

where  the  non-Bessel  factor  0(/3)  is  finite  in  the  range  of  integration,  and  asymptotically 

m ' ar"  (3.2) 

in  which  a is  constant.  The  value  of  n in  our  formulas  in  the  previous  section  is  = 4.  The  general 
case  is  discussed  in  Appendix  I. 

We  will  replace  the  infinite  integral  (3.1 ) with  a finite  integral.  Given  a large  value  N,  we  can  esti- 
mate an  upper  bound  of  the  absolute  Integral, 

f \m  f^mJom\dp  0.3) 

N 
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called  the  absolute  remainder,  by  substituting  the  asymptotic  expansions  of  ipifi)  and  Bessel  functions. 
We  let  the  trigonometric  functions  in  the  latter  equal  one.  Denoting  the  absolute  remainder  by 
[Abs  /)  N. 

[4bs/|N  < [2(7/(7rv^)l  (nA'")-’  . (3.4) 


Let  e be  the  error  we  can  tolerate  in  our  computation.  In  our  actual  computation,  we  chose 
e = 10-5. 

The  value  of  N is  evaluated  by  equating  the  right  hand  side  of  (3.4)  to  e; 

[2o/(7r\/oT)|  (n/V")-’  = e.  (3.5) 

Then,  integral  / in  (3.1)  is  approximated  by 

N 

/ =.  / m y,(/i4)  Jo(^R)  d^-  (3-6) 

0 

The  value  of  A'  was  small  in  most  of  our  computation  | A = 10  except  in  (6.4)  ] , and  our  numerical 
scheme  worked  very  effectively. 

We  list  in  the  following  the  asymptotic  expansions  of  the  non-Bessel  factors  0(^)  contained  in  the 
integrands  of  our  integral  solutions  (2.16),  (2.19),  and  (2.20); 

a,  --  IT'* 

02  ~ t(1  +(r'’) 

e'-iT  ~i 


(T-a, )/v/5£SC  ~ /r^-^8 

(T-a2)/VDESC  ~ -(1  -£)r®. 


4.  RAMP/STEADY  LOADING 

We  used  two  load  tests  to  fit  our  theoretical  curves.  One  was  the  Sun  Oil  Corporation’s 
(SUNOCO)  data  obtained  during  the  winters  of  1973-1974  and  1974-1975  at  Resolute  Bay,  North- 
west Territory  (unpublished).  The  other  was  Frankenstein’s  data^  obtained  on  Portage  Lake, 
Michigan,  and  the  Garrison  Dam  Reservoir,  North  Dakota,  on  20  March  1956  and  18  January  1957. 
Among  these  tests,  we  chose  the  ideal  ramp/steady  loading  for  our  numerical  computation.  In  this 
loading,  as  illustrated  in  Figure  2,  the  load  P is  increased  initially  at  a constant  rate  P and,  after  a 
certain  time  Tq,  kept  constant  at  £ = PTq.  However,  since  SUNOCO  does  not  allow  the  publication 
of  their  data,  we  cannot  include  their  data  in  this  paper. 

We  will  derive  the  ramp/steady  formulas  by  use  of  the  step-loading  formulas  given  in  the  previous 
section.  However,  since  both  SUNOCO  and  Frankenstein  measured  only  deflection,  we  derive  only 
the  deflection  formulas. 

Define  the  influence  function  wq(T)  by  letting/’  = 1 in  (2.16): 
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Substituting  (4.1 ) into  (4.2)  and  integrating  with  regard  to  X,  we  get  the  deflection  w(T)  for 


w[T)  = [PKiiApQ})]  (fy,  - (^2  + .73)  (4.4) 

where 

N 

Uy  = f -L  (e-"’^-! +a,  n /i(^)  /o(^^)  (4-5) 

0 “1 

73  = / (1  -e'“2T)  y,{^n)  y (^/?)  (4.7) 

^ y/DESC  «2 

The  absolute  remainders  are  as  follows: 

Mftsfy,]”  < [^^|{4n^/AR)]  rr*  (4.8) 

[/16s  6^2 1“  < [T/{2ns/M)]N-‘*  (4.9) 
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i/tfestyjj;;;  < [(i  -m2K^)]  J-  (i  . (4.10) 

Substituting  (4.1)  into  (4.3)  and  integrating  with  regard  to  X,  we  get  the  deflection  w{T)  for 
To  i T : 

w(7)  = (/’/(7r.4p«2)](/, +/2+/3)  (4.11) 

where 

/,  = / 1 - ^ (e“t  ^0  _ 1 ) ) /g(0«)  d&  (4.1 2) 

/2  = / (1  -e-“l^)^  (e“>^0-1)yi(/i4)  J^m)  d&  (4.13) 

/3  = / e-“2T  ^ (e“2To_i)  y^((i4)  ]^{fiR)  d&.  (4.14) 


The  absolute  remainders  are  as  follows; 

(/16s/,  I “ < (2;rv/4^)'’ 

[/16s /jl”  < [TI(2rs/M)]N-* 

(/I6s/3|“  < |(?-f)/(2frv^)){(e^^O_])/(r7-o)]  e'^'^N-^. 
Computer  programs  for  these  formulas  are  shown  in  Appendix  II. 


(4.15) 

(4.16) 

(4.17) 


5.  CURVE  FITTING  TO  TIME  LAPSE  DEFLECTIONS 

Frankenstein^  placed  a 1 2-ft-diameter  tank  on  the  ice  and  pumped  the  adjoining  water  into  the 
tank.  (We  call  this  the  distributed  load  test.)  However,  the  temperature  of  the  water  in  the  tank 
obviously  disturbed  the  ice  temperature.  He  then  tried  a variation  by  placing  a 1 7.3-in.-diameter 
concrete  block  under  the  1 2-ft-diameter  tank.  (We  call  this  the  concentrated  load  test.).  The  water 
in  the  tank  was,  in  this  case,  isolated  from  the  ice  and  did  not  disturb  the  ice  temperature. 

The  load-vs-time  curves  of  these  tests  and  the  measured  deflections  are  shown  in  Figures  3 and 
4.  "TANK”  designates  the  deflection  of  the  edge  of  the  tank.  “RODS”  are  the  sites  where  the 
measurements  were  taken.  The  distances  of  the  measurement  sites  from  the  center  of  the  load  are 
listed  in  Tables  I and  II. 

The  material  constants  found  by  the  curve  fitting  are  shown  in  Tables  I and  II.  They  vary  with 
the  location  of  the  measurement. 

To  show  the  significance  of  the  material  constant  variation  with  the  measurement  sites,  we  chose 
the  material  constants  determined  at  rod  1 of  Frankenstein’s  concentrated-load  time-lapse  curve,  and 
computed  the  deflections  at  the  other  measurement  sites.  Figure  5 shows  the  comparison  of  the 
computed  curves  and  the  measured  data.  The  left  and  right  columns  show  the  ramp  and  steady 
portions  of  the  deflection  curves,  respectively.  They  are  designated  by  (r)  and  (s)  respectively 

To  express  the  degree  of  curve  fitting  we  devised  the  trapezoidal  error  (TE).  In  Figure  6,  a,  A 
and  b,  B show  two  pairs  of  measured  and  computed  deflections  at  two  consecutive  times  /,  and  /j. 
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Figure  4.  Concentrated  load  test  by  Frankenstein  (ref.  3,  test  8). 


Distonce  I m) 


Figure  6.  Elements  of  TE.  Figure  7.  The  TE  of  Frankenstein 's  distributed  load  test 

(ref  3,  test  5). 


respectively.  We  squared  A-a  and  B-b  and,  in  case  of  the  upper  figure  where  the  errors  are  of  the 
same  sign,  computed  the  area  of  the  trapezoid  of  which  the  bases  are  (A-a)^  and  (B-b)^  and  the 
height  ^2  “ case  of  the  lower  figure  where  the  errors  change  sign,  we  calculated  the  sum  of 
the  areas  of  the  two  triangles  of  which  the  bases  are  (A-a)^  and  {b-B)^  and  the  heights  /q  - ^ j and 
^2  " ^0  tespectively,  where  tQ  is  the  abscissa  of  the  intersection.  Denoting  by  5 the  area  of  such  a 
figure,  we  defined  TE  by 

TE=V(2W  (5.1) 

where  the  summation  is  over  all  the  intervals  and  T the  sum  of  the  abscissa  intervals. 

The  TE  indicates  a sort  of  absolute  maximum  error.  Its  unit  is  m.  If  the  deflections  are  of 
ordinary  magnitude,  the  TE  of  order  lO'^  and  1 0'^  means  a good  and  tolerable  fit,  respectively.  If 
the  deflections  are  very  small,  as  in  the  case  of  rod  4,  the  smallness  of  the  value  of  TE  docs  not  mean 
much.  We  did  not  list  the  computed  values  at  rod  4 in  Tables  I and  II. 

We  evaluated  the  TE  for  all  the  possible  cases.  They  are  shown  in  Figures  7 and  8.  The  abscissa 
is  the  distance  from  the  center  of  the  load.  The  measurement  sites  are  noted  on  the  abscissa  axis. 

The  circled  points  are  those  whose  material  constants  are  used  to  compute  a set  of  TE.  The  sets  of 
TE  thus  computed  are  connected  with  solid  or  broken  lines  and  labeled  with  the  appellations  of  the 
circled  measurement  sites. 

Comparison  of  Figures  7 and  8 shows  that  the  concentrated  load  test  has  smaller  overall  TE  values 
than  the  distributed  load  test.  Flowever,  we  cannot  recognize  any  significant  effect  of  the  tempera- 
ture distribution  due  to  the  watertank  temperature  disturbance.  Probably  the  cracks,  whose  appear- 
ances are  noted  in  Figures  3 and  4 but  are  not  considered  in  our  formulation,  were  more  detrimental 

After  the  numerical  computation  was  finished  in  1976,  Dr.  Andrew  Assur,  an  ice  mechanics 
expert  at  CRREL,  notified  us  that  the  variations  of  material  constants  in  Tables  I and  II  are  in  the 
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range  of  reasonable  values  from  the  viewpoint  of  the  nonlinear  viscoelastic  constants  (Shumskij'^ ). 
We  tried  in  1977  to  reevaluate  the  material  constants;  we  thought  that,  although  the  theoretical 
curve  is  formulated  on  the  linear  assumption,  if  we  fit  the  theoretical  curve  in  the  narrow  time 
interval  and  space  span,  we  can  find  the  material  constants  close  to  the  incremental  viscoelastic 
constants.  However,  this  plan  could  not  be  executed  because  the  distances  between  the  measuring 
rods  were  too  large. 


6.  ASYMPTOTIC  DEFLECTION 


We  shall  show  in  the  following  that  only  one  material  constant  is  contained  in  the  asymptotic 
formulas.  The  curve  fitting,  therefore,  must  be  carried  out  in  the  initial  stage. 

Referring  to  the  asymptotic  relationships  in  (3.7),  we  find  that,  when  T is  large,  both  the  step- 
loading formulation  (2.16)  and  the  ramp/steady  loading  formulation  (4.1 1 ) reduce  to 

w = [/>/(rr/lp82)l  (K//i)  (6.1) 

where 

oo 

= J (1 -e-T0‘^)y,  (M)yo(l3«)^'/5-  (6.2) 

0 

It  is  assumed  in  this  derivation  that  t > 0,  and  that  only  large  values  of  (i  are  effective  in  the  integra- 
tion. Letting 


X = pA 

(6.2)  becomes 


(6.3) 


K=  f (\-e^^^'^)J^{x)Jol(R/A)x]dx 
0 


(6.4) 


where 

Ta  = TA^ 

= tat  12/a(1  -v) 

n\ 


(6.5) 

(6.6) 


Thus,  all  the  material  constants  are  lumped  into  the  second  factor  of  (6.6).  The  stress  formulas,  al- 
though not  mentioned  here,  can  be  similarly  transformed. 

As  shown  in  Appendix  I,  (6.4)  cannot  be  analytically  integrated;  it  must  be  numerically  integrated. 
To  effect  the  numerical  integration,  the  non-Bessel  factor  in  (6.4)  is  so  chosen  that  it  becomes  zero 
atx  = The  absolute  remainder  is  estimated: 


Mbs /Cl”  < 


(6.7) 


Graphs  of  integral  K for  the  values  of  R/A  = 0.2  and  2.0  are  shown  in  Figure  9.  When  Ta  = “, 
the  non-Bessel  factor  becomes  equal  to  one.  At  this  limit,  therefore,  K = 1 when  R < /4,  and  K = 

0 when  R > A.  As  shown  in  the  graphs,  this  limit  is  almost  reached  when  Tf^  > 1000. 

Exact  integral  K was  formulated  for  the  ramp/steady  loading,  and  evaluated  by  use  of  a set  of 
consunts:  Tq  - 6x10^  sec,  r = 10,  £ = 1/6,  T?i/fo  = 6.12x10^  sec  = 17hr,  £(,  = 7x10® 
kg/m2,r-  = 0.5,and/4  = 0.5.  These  consUnts  give  8 = 29.31  m and  Ta  = f(2.48x10''’ day'). 

As  shown  in  Figure  9,  the  asymptotic  integral  K is  very  close  to  the  exact  integral  in  the  range 

Ta  > 0.1 . The  above  constants  are  the  rough  estimates  used  before  starting  the  elaborate  calculations. 
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They  are  not  listed  in  the  Tables.  We  did  not  use  other  sets  of  constants  to  evaluate  the  exact 
integral  K.  We  expect  that  all  the  exact  curves  should  show  the  similar  coincidence  with  the 
asymptotic  curve  although  with  individual  variations. 

The  values  of  at  the  final  time  of  the  two  tests  are  listed  in  Table  III.  These  values  are  very 
small.  However  we  experienced  that  the  modification  of  some  material  constants  was  insensitive  on 
the  modification  of  the  computed  deflection  values. 


Table  III.  Final  time  of  the  three  tests. 


In  physical  unit 

In  T/i,  unit 

Frankenst'cn's  distributed  ioad  test 

420  min  = 7 hr 

2.67X10'^ 

concentrated  load  test 

240  min  = 4 hr 

1.2XI0‘* 

7.  DEFLECTION  PROFILES 

We  computed  the  deflection  profiles  of  the  concentrated  load  test  (Frankenstein,^  test  8)  at  12.7, 
32.6,  and  11 8.5  min  by  use  of  the  material  constants,  T = 10,  E = 0.02,  = 7x10^  sec,  and 

Eq  = 2x10®  kg/m^,  as  shown  in  Figure  10.  These  material  constants  are  round  numbers  intermediate 
between  the  material  constants  at  rod  1 and  rod  2 in  Table  II.  The  three  chosen  times  mentioned 
above  are  marked  in  Figure  4.  The  computed  profiles  are  quite  different  from  the  measured  profiles. 
We  varied  the  material  constants  but  could  not  find  values  that  make  the  theoretical  curve  assume  a 
similarity  to  the  measured  curve.  It  is  our  impression  that  the  measured  profiles  do  not  belong  to 
the  family  of  curves  that  our  analytical  formula  can  describe.  The  measured  and  computed  curves 
intersect  between  rod  1 and  rod  2,  indicating  the  reliability  of  our  computation,  as  may  be  expected 
from  the  choice  of  the  material  constants. 


Figure  10.  Deflection  profile.  (Frankenstein's  concentrated  load  test,  ref.  3, 
test  8.) 
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Figure  11.  Asymptotic  deflection  profile.  (Theoretical  continuation 
of  Frankenstein's  concentrated  had  test,  ref.  3,  test  8.) 


At  f = the  integral  K in  (6.4)  becomes 
K = 1 for  0 < /•  < (7 

(7.1) 

= 0 for  0 < a < /•. 

The  deflection  at  / = is 

= q/p  for  0 < r < a 

(7.2) 

= 0 for  0 < a < /•. 

Therefore,  the  water  tank  sinks  theoretically  to  w„  = 93.3  m in  the  case  of  the  concentrated  load. 
However,  the  ice  thickness  h is  0.556  m.  Our  analytical  formulas,  therefore,  become  invalid  beyond 
a certain  elapsed  time.  (In  the  case  of  the  distributed  load  test,  w„  = 1.350  m and/)  = 0.597  m.) 

Theoretical  deflection  profiles  for  large  times  are  shown  in  Figure  11 . At  time  infinity,  our 
analytical  deflection  comes  to  the  vertical  line  denoted  by  r = «>.  Because  5.41  x1  O'®  xr  (days) 
in  the  case  of  the  concentrated  load  test,  the  largest  time,  500  days,  chosen  for  this  calculation  is  still 
voo  short.  However,  the  mode  of  approach  to  the  ultimate  t = °°  curve  is  observable  with  the  curves 
mFigurell.  (In  the  case  of  the  distributed  load  test,  = 1.98x10*^xf  (days).] 
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APPENDIX  I.  ANALYTICAL  BACKGROUND 


% 

1 


A.  The  following  theorem  shows  the  condition  under  which  the  integral  (3.1)  becomes  either 
discontinuous  or  continuous  at /?  = A. 

Theorem  1.  The  integral  (3.1)  is  discontinuous  or  continuous  HR  = A when  n in  (3.2)  is  equal 
to  or  larger  than  zero,  respectively. 

Proof.  We  can  rewrite  (3.1 ) to  a one-parameter  integral 

/(a)  = f(x,  a)  dx  (A.1 ) 

0 

bylettingx  = /i4,i.e.  a = where  f(x,  a)  is  continuous  with  regard  to x and  a.  The  condition 
that  /(a)  is  a continuous  function  of  a is  that  the  integral  (A.1 ) converges  uniformly  with  respect  to 
a (c.f.  Titchmarch,'®  p.  25).  The  integral  (A.1)  uniformly  converges  when  n > 0,  but  does  not 
when  /7  = 0. 

B.  We  shall  consider  in  the  following  the  integral  (3.1 ) whose  non-Bessel  factor  0(/3)  is  finite  in  the 
range  of  integration  but  asymptotically  becomes  zero  in  a more  general  form  than  in  the  specific 
form  (3.2). 

Let  an  asymptotic  expansion  of  (pifi)  be 

«(^)  ' ^ 0n((5).  (B.1) 

Rewrite  (3.1 ) as 

/ = /o  + ^ (B.2) 

where 


^0 

and 


/,(/3-4)  jom 


(B.3) 


] 


We  choose  such  an  integer  m that  makes  /q  rapidly  convergent.  We  choose  such  a function 
that  makes  (B.4)  analytically  integrable.  The  following  theorem  is  useful  for  the  choice  of 

Theorem  2.  Let  F(z)  be  an  even  function  of  the  complex  variable  z = x + iy  that  becomes  zero 
atz  = “>  and  possesses  only  algebraic  singularities  (pole  or  branch  points)  on  the  upper  half  plane 
but  no  poles  on  the  real  axis.  Then 

oo 

j F(x)  y,(</x)  ]Q[bx)  dx 

0 

+ ep 

= 5^^)  W^’)(az)  jQ(bz)dz  when  r?  > b > 0 

• \jO 

(B.S) 

+ OP 

= ^ ^ F(z)  J , (az)  ) {bz)  dz  when  Q < a < b 

— PP 

(B.6) 

+ OP 

where  dz  means  the  integral  along  the  contour  in  Figure  12,  where  radius  e is  infinitesimal,  and 

the  z-plane  is  cut  along  the  negative  real  axis. 

Proof.  Consider  the  contour  integrals 

4-  oo 

/(a  > b)  = ^ Kjj  F(z)  (az)  /o(bz)  dz 

(i) 

where  a > b > 0,  and 

l{a  < b)  = ^ ^ F{z)  ;,(az)  ^(j^Hdz)  dz 

(ii) 

when  0 < a < b.  Use  of  the  asymptotic  formulas  show  that  (az)  Jq 

(bz)  are  zero  on  the  infinitely  large  circle  when  a > b > OandO  < a 
Therefore  we  may  consider  only  the  contour  along  the  real  axis 

(bz)  and  ) , (az) 

< b,  respectively. 

4 oo 

/(a  > b)  = i y Flz)  > (az)  /g  (bz)  dz 

(iii) 

4 oo 

/(a  < b)  = 1 J F(z)  y,  (az)  (bz)  dz. 

(iv) 
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We  divide  the  real  axis  in  three  parts; 
region  - ~ - e,  and  z = x in  the  region  e 


e,-e  ~ e,  and  e ~ We  let  z = - x in  the 

~ neglect  the  infinitesimal  terms,  and  let 


F(z)  = F(0) 

) (az)  = - 2J/(naz) 

/yy)(/>z)  = [2/7(6z)llog(6z/2). 

Then  (iii)  and  (iv)  become 

oe 

/(a  > b)  = -iF(0)+  f F(x)  y,(ax)  Jo{bx)  dx  (v) 

0 

oo 

/(a  < b)  = f F(x)  y,(ax)  Jo{bx)  dx.  (vi) 

0 

Equations  (v)  and  (vi)  prove  (B.5)  and  (B.6),  respectively. 

C.  The  need  of  Theorem  2 appears  frequently  in  the  mathematical  study  of  the  problems  of  a float- 
ing ice  plate  and  the  problems  of  an  elastic  plate  on  an  elastic  foundation.  A similar  integral  in- 
cluding only  one  Bessel  function  was  proved  by  Dougal  (ref.  1,  p.  138  and  147)  as  early  as  in  1903. 

When  f = 0,  our  solution  of  the  viscoelastic  plate  reduces  to  the  solution  of  the  elastic  plate. 

The  elastic  solution  thus  found  is  composed  of  the  following  integrals; 


M 


oo 

0 = f — ~ h (<^^)  hibx)  dx 

/- 

K 1 +. 


= / 


1 +x‘' 


y,  (C7X)  y,  (bx)  dx 


J]  (ox)  Joibx)  dx 


where  (7  = AF^andi)  = FF^.  We  can  carry  out  these  integrals  by  direct  or  indirect  application 
of  Theorem  2; 

Mn 


berb  ker'o  - beib  kei'o  + a'^ 

when  b ^ 

ber'a  kerb  - bei'o  keib 

when  a = 

- ber'b  ker  o + bei'b  kei'o 

when  b ^ 

- ber’a  ker'b  + bei’o  kei'b 

when  a 5 

beib  ker'o  + berb  kei'a 

when  b ^ 

bei'o  kerb  + ber’a  keib 

when  a ^ 
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Mq  and  Afj  are  found  by  directly  applying  the  theorem.  is  found  by  differentiating  Mq  with 
regard  to  B.  Wyman”  derived  Mq  by  integrating  a concentrated-load  elastic-plate  solution  over  the 
loading  circle. 

The  continuity  of  Mq  and  = b is  obvious  on  the  strength  of  Theorem  1 . We  shall  show, 

however,  a direct  proof  in  the  following.  We  shall  prove  that 

berA’  kerV  - beiA  keiV  + at*’  = ber'x  kerx  - bei'A  kelx  (C.1 ) 


and 


beix  ker'x  + berx  kei'x  = bei'x  kerx  + ber'x  keix. 
To  prove  this,  note  that 

w,  (x)  = berx  + / beix 

and 

tvj  (x)  = kerx  + / keix 
are  the  solutions  of  the  differential  equation 

= 0. 

^^2  X ax 

This  can  be  proved  by  decomposing  the  equation 

Wx^  X o-x/ 

of  which  (C.3)  and  (C.4)  are  the  solutions. 

We  can  find  that  the  Wronskian 

w,  (x)  W2(x) 

/,  (x)  W2(x) 

is  equal  to 

= 


Thus  we  have  the  identity 

berx  + /'  beix  kerx  + /'  keix 
ber'x  + / bei'x  ker'x  + / kei'x 


= _i 

X 


of  which  the  real  part  gives  (C.1 ) and  the  imaginary  part  gives  (C.2). 
Theorem  2 can  be  extended  in  many  ways.  Nevel  '*  found  that 


» -f  «e 

J F{x)  dx  = ^ f F(z)  logz  dz 


(C.2) 


(C.3) 


(C.4) 


(C.3) 
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for  an  odd  function  F{z)  that  does  not  have  any  pole  on  the  real  axis  and  vanishes  atz  = <*>. 

D.  It  is  impossible  to  apply  Theorem  2 to  the  integrals  of  w in  (2.16),  in  (2.19),  and  Og  in  (2.20) 
for  the  following  reason. 

The  function  exp(-  aj  0 has  essential  singularities  at  the  roots  of 


^ + £ = 0, 

because 


lim  02  = «> 
0^-»-£. 


The  function  exp(-o,  T)  does  not  possess  any  essential  singularities  because  the  limit  of 
Q,  = 2t/[t^  + 1 + t + y/(T^  + 1 + T)i  - 4t{^  + f)  ] 


is  finite.  However,  the  real  part  of  a,  becomes  negative,  and  exp(-a,  T)  diverges,  as  | ^ | -*■»<>  in  a 
certain  range  of  direction. 

Theorem  2 does  not  apply  to  integral  K in  (6.4)  because  the  point  x =0  is  an  essential  singularity. 
The  only  alternative  we  can  find  for  the  integration  of  (3.1)  is  the  use  of  Barnes’  integral  method. 

It  consists  in  substituting  the  integrals 


/ M = J_  ?'  r(-5)(’/2x)v^2S 
’ 2m  J,  r(v  +S  +1)  " 


(D.1) 


Ml)  (z)  = -L 

''  2iii 


-c+«>i  / 

/ r(-v-s)r(-s)(-i 

'Oej  \ 


z\  ds 


(D.2) 


for  l^(x)  and  respectively,  where  c is  a real  number  satisfying  c > R(v),z  is  complex,  and 

X is  real.  We  can  usually  exchange  the  order  of  integration  to  carry  out  the  integration  with  regard 
to  X or  z.  Then,  we  can  carry  out  the  rest  of  the  integration  in  most  cases  by  use  of  the  theorem 
of  residue.  Only  the  forms  (D.1 ) and  (D.2)  serve  this  purpose.  The  other  Barnss’  representations 
of  (x ) and  ) (z)  do  not  enable  us  to  carry  out  the  above  two  procedures. 

However,  as  mentioned  by  Watson  (ref.  1 8,  p.  1 92),  (D.1 ) does  not  hold  true  for  v = 0,  and 
(D.2)  does  not  hold  true  when  v = 0 and  z is  real.  In  these  two  cases,  the  integrands  of  (D.1 ) and 
(D.2)  become  proportional  to  s'l  as  s approaches /o®  as  the  limit  on  the  imaginary  axis.  Therefore 
we  cannot  use  Barnes’  integral  method  to  carry  out  our  integrals. 


APPENDIX  II.  COMPUTER  PROGRAMS 
Ramp  Time  Profiles 
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